Show the code
pacman::p_load(sf, tidyverse, questionr, janitor, psych, ggplot2, gcookbook, tmap, ggpubr, egg, corrplot, gtsummary, regclass, caret, heatmaply, ggdendro, cluster, factoextra, spdep, ClustGeo, GGally, skimr, stringr)Regionalisation with Spatially Constrained Cluster Analysis
case study : Regionalisation by Multivariate Water Point Attributes with Non-spatially Constrained and Spatially Constrained Clustering Methods.
1. OVERVIEW
Regionalisation with Spatially Constrained clustering analysis requires similar observations to be grouped according to their statistical attributes and spatial location.
1.1 Objectives
1.2 Study Area Introduction
1.3 Scope of Works
2. R PACKAGE REQUIRED
The following are the packages required for this exercise :
2.1 Load R Packages
3. GEOSPATIAL DATA
3.1 Acquire Data Source
Aspatial Data
- Download the Nigeria data set in shapefile format via Access WPdx+ Global Data Repository from WPdx Global Data Repositories.
- Rename the title of the data set to “geo_export”.
Note
The file size of the downloaded data is about 422 MB due to water points data from multiple countries.
- Such file size may require extra effort and time to manage the code chunks and files in the R environment before pushing them to GitHub.
Hence, to avoid any error in pushing files larger than 100 MB to Git, filtered Nigeria water points and removed unnecessary variables before uploading into the R environment.
Therewith, the CSV file size should be lesser than 100 MB.
Geospatial Data
- Download the Nigeria geoBoundaries data set at ADM2 level4 from geoBoundaries.org or the Humanitarian Data Exchange portal.
- Rename the title of the data set to “nga_admbnda_adm2_osgof_20190417”
4 Runfola, D. et al. (2020) geoBoundaries: A global database of political administrative boundaries. PLoS ONE 15(4): e0231866. https://doi.org/10.1371/journal.pone.0231866
3.2 Import Attribute Data
Four (4) data frames to be created for different context, i.e.
wp_coord = coordinated related variables.
wp_cond = status and conditions related variables.
wp_adm = administrative and measurements related variables.
wp = master file that combine all three (3) data frames above.
3.2.4 Create Master File
3.2.5 Convert Well Known Text (WKT) Data to SF Data Frame
The “New Georeferenced Column” in wp_rds contains spatial data in a WKT format.
Two (2) steps to convert the WKT data format into an sf data frame.
3.2.5.1 derive new field :: “geometry”
3.2.5.2 convert to SF Data Frame
3.3 Import Boundary Data
::: {.callout-warning appearance=“simple” icon=“false”} ::: {.callout-alert appearance=“simple” icon=“false”}
Usage of the code chunk below :
st_read( ) - sf - to read simple features.
select( ) - dplyr - to select “shapeName” variable. :::
Show the code
bdy_nga <- st_read(dsn = "data/geospatial",
layer = "geoBoundaries-NGA-ADM2",
crs = 4326) %>%
select(shapeName)
problems(bdy_nga)3.3.1 save and read RDS file :: wp_adm
3.3.2 Review Imported File
3.3.2.1 check for missing and duplicated data
3.3.2.2 list the unique “shapeName” associated with duplication
3.3.2.3 verify findings in section 3.3.1.2
3.3.1.4 acquire State info for duplicated value
| lga | row_id | headquarter | state | iso3166code | dupl_shapeName_coord | lga_office_coord |
|---|---|---|---|---|---|---|
| Bassa | 94 | Oguma | Kogi | NG.KO.BA | 7.791971, 7.031827 | 7.80932, 6.74853 |
| Bassa | 95 | Bassa | Plateau | NG.PL.BA | 10.08015, 8.782946 | 10.11143, 8.71559 |
| Ifelodun | 304 | Share | Kwara | NG.KW.IF | 8.544311, 5.052235 | 8.5 5.0 |
| Ifelodun | 305 | Ikirun | Osun | NG.OS.ID | 7.920948, 4.636735 | 7.5 4.5 |
| Irepodun | 355 | Omu Aran | Kwara | NG.KW.IR | 8.169349, 4.926215 | 8.5 5.0 |
| Irepodun | 356 | Ilobu | Osun | NG.OS.IP | 7.84861, 4.498797 | 7.5 4.5 |
| Nasarawa | 519 | Bompai | Kano | NG.KN.NA | 12.00446, | |
| 8.578262 | 11.5 8.5 | |||||
| Nasarawa | 520 | Nasarawa | Nasarawa | NG.NA.NA | 8.304034, 7.760272 | 8.53477, 7.70153 |
| Obi | 546 | Obarike-Ito | Benue | NG.BE.OB | 7.022495, 8.281026 | 7.01129, 8.33118 |
| Obi | 547 | Obi | Nasarawa | NG.NA.OB | 8.35534, 8.734777 | 8.37944, 8.78561 |
| Surelere | 693 | Surelere | Lagos | NG.LA.SU | 6.493217, 3.346919 | 6.50048, 3.35488 |
| Surelere | 694 | Iresa-Adu | Oyo | NG.OY.SU | 8.088897, 4.393574 | 8.08459, 4.38538 |
3.4 Data Wrangling
3.4.1 Edit Duplicated Value :: “shapeName”
3.4.1.1 validate edited value :: “shapeName”
3.4.2 Perform Point-in-Polygon Overlay
Combine both attribute and boundary of the water points into a simple feature object.
3.4.2.1 join objects :: wp_sf, bdy_nga
3.4.2.2 inspect joined file :: wp_joined
-- retrieve geometry summary
-- assess uniqueness of Water Point
-- determine reference point :: “shapeName” or “clean_adm2”
3.4.3 Remove Water Point Outside LGA Boundary
3.4.4 Replace “NA” with “Unknown”
3.4.4.1 identify variable with missing value
3.4.4.2 replace “NA” with “Unknown”
3.4.5 Combine Unique Value
3.4.6 Extract Functional Water Point for New Table :: wpt_functional
3.4.6.1 inspect variable and value
-- reveal value :: “crucialness_score”
-- water points over usage limit
3.4.6.3 Exploratory Data Analysis (EDA) :: wpt_functional
-- plot “status_clean”
Show the code
ggplot(data = wpt_functional,
aes(fct_infreq(status_clean), fill=status_clean))+
geom_bar()+
geom_text(
aes(label=after_stat(count)),
stat='count',
nudge_x=-0.25,
vjust=-0.2)+
geom_text(
aes(label= scales::percent(signif(after_stat(count/sum(count))))),
stat='count',
nudge_x=0.25,
vjust=-0.2)+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
guides(fill=guide_legend (title = 'Status'))
-- plot “water_tech_category”
Show the code
ggplot(data=wpt_functional,
aes(x=fct_infreq(
water_tech_category)))+
geom_bar(aes(
fill = water_tech_category),
width = 0.8)+
geom_text(aes(
label = ..count..),
stat = "count",
vjust=-0.2,
size = 3.5,
color = "black")+
scale_x_discrete(guide = guide_axis(n.dodge = 2))+
guides(fill=guide_legend (title = 'Water Tech Deployed'))Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
ℹ Please use `after_stat(count)` instead.
-- plot “water_source_clean”
Show the code
ggplot(data=wpt_functional,
aes(x=fct_infreq(
water_source_clean)))+
geom_bar(aes(
fill = water_source_clean),
width = 0.8)+
geom_text(aes(
label = ..count..),
stat = "count",
vjust=-0.2,
size = 3.5,
color = "black")+
scale_x_discrete(guide = guide_axis(
n.dodge = 2))+
guides(fill=guide_legend (
title = 'Source of Water Supply'))
3.4.5.4 add attribute to new data table
3.4.5.5 extract non-functional water point
3.4.5.6 inspect variable and value
3.4.5.7 EDA :: wpt_nonFunctional
-- plot “status_clean”
Show the code
ggplot(data = wpt_nonFunctional,
aes(fct_infreq(status_clean),
fill=status_clean))+
geom_bar()+
geom_text(
aes(label=after_stat(count)),
stat='count',
nudge_x=-0.25,
vjust=-0.2)+
geom_text(
aes(label= scales::percent(
signif(
after_stat(
count/sum(count)
)))),
stat='count',
nudge_x=0.25,
vjust=-0.2)+
scale_x_discrete(
guide = guide_axis(
n.dodge = 2))+
guides(fill=guide_legend (
title = 'Status'))
-- plot “water_tech_category”
Show the code
ggplot(data=wpt_nonFunctional,
aes(fct_infreq(
water_tech_category)))+
geom_bar(aes(
fill = water_tech_category),
width = 0.8)+
geom_text(aes(
label = ..count..),
stat = "count",
vjust=-0.2,
size = 3.5,
color = "black")+
scale_x_discrete(guide = guide_axis(
n.dodge = 2))+
guides(fill=guide_legend (
title = 'Water Tech Deployed'))
-- plot “water_source_clean”
Show the code
ggplot(data=wpt_nonFunctional,
aes(fct_infreq(
water_source_clean)))+
geom_bar(aes(
fill = water_source_clean),
width = 0.8)+
geom_text(aes(
label = ..count..),
stat = "count",
vjust=-0.2,
size = 3.5,
color = "black")+
scale_x_discrete(guide = guide_axis(
n.dodge = 2))+
guides(fill=guide_legend (
title = 'Source of Water Supply'))
3.4.5.8 add wpt_nonFunctional to wp_nga
3.4.5.9 extract unknown water point
3.4.5.10 inspect variable and value
-- reveal value :: “status_clean”
3.4.5.11 add wpt_unknown to wp_nga
3.4.5.12 visualise distribution :: “status_clean”
Usage of the code chunk below :
qtm( ) - tmap - to plot a thematic map quickly.
tmap_arrange( ) - tmap - to arrange small multiples in grid layout.
Show the code
total_wp <- qtm(wp_nga,"total_wp")+
tm_layout(legend.height = 0.3, legend.width = 0.5)
wp_functional <- qtm(wp_nga,"wp_functional")+
tm_layout(legend.height = 0.3, legend.width = 0.5)
wp_nonFunctional <- qtm(wp_nga,"wp_nonFunctional")+
tm_layout(legend.height = 0.3, legend.width = 0.5)
wp_unknown <- qtm(wp_nga,"wp_unknown")+
tm_layout(legend.height = 0.3, legend.width = 0.5)
tmap_arrange(total_wp, wp_functional, wp_nonFunctional, wp_unknown, asp=0, ncol = 2, nrow = 2, widths = 5, heights = 10, sync = TRUE)
3.4.5.13 extract “water_tech_category” to wp_nga
3.4.5.14 visualise wp_nga distribution :: “water_tech_category”
Show the code
handPump <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_handPump",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
mechPump <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_mechPump",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
tapStand <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_tapStand",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
tmap_arrange(handPump, mechPump, tapStand,
asp=1,
ncol=2,
sync = TRUE)
3.4.5.15 extract “usage_capacity” to wp_nga
3.4.5.16 visualise wp_nga distribution :: “usage_capacity”
Show the code
uc300 <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_uc300",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
uc1000 <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_uc1000",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
uc250 <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_uc250",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
tmap_arrange(uc300, uc1000, uc250,
asp=1,
ncol=2,
sync = TRUE)
3.4.5.17 extract “is_urban” to wp_nga
3.4.5.18 visualise wp_nga distribution :: “is_urban”
Show the code
urban1 <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_urban1",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
urban0 <- tm_shape(bdy_nga)+
tm_polygons(alpha = 0.1) +
tm_shape(wp_nga) +
tm_dots(col = "pct_urban0",
border.col = "gray60",
border.lwd = 0.5) +
tm_view(set.zoom.limits = c(5,9))
tmap_arrange(urban1, urban0,
asp=1,
ncol=2,
sync = TRUE)
3.4.5.19 save and read RDS File :: wp_nga
Show the code
write_rds(wp_nga,"data/geodata/wp_nga.rds")
wp_nga <- read_rds("data/geodata/wp_nga.rds")3.4.5.20 transform to Projected Coordinate System
3.5 Exploratory Data Analysis
3.5.1 Identify Outliers
3.5.1.1 plot boxplot “pct_functional”
Show the code
ggplot(data=wp_ngaTrans,
aes(x=`pct_functional`)) +
geom_boxplot(color="black",
fill="#543005")
3.5.1.2 plot boxplot “pct_nonFunctional”
Show the code
ggplot(data=wp_ngaTrans,
aes(x=`pct_nonFunctional`)) +
geom_boxplot(color="black",
fill="#C16622FF")
3.5.1.3 plot boxplot “pct_unknown”
Show the code
ggplot(data=wp_ngaTrans,
aes(x=`pct_unknown`)) +
geom_boxplot(color="black",
fill="#FFA319FF")
Remarks :
Among these 3 key categories of “status_clean”, “unknown” has the most outliers.
3.5.2 Multi-plot Histogram
3.5.2.1 plot histogram for “status_clean”
Show the code
pctFunctional <- ggplot(data = wp_ngaTrans,
aes(x = `pct_functional`))+
geom_histogram(bins=10,
colour = "black",
fill = "#543005")
pctNonFunctional <- ggplot(data = wp_ngaTrans,
aes(x = `pct_nonFunctional`))+
geom_histogram(bins=10,
colour = "black",
fill = "#C16622FF")
pctUnknown <- ggplot(data = wp_ngaTrans,
aes(x = `pct_unknown`))+
geom_histogram(bins = 10,
colour = "black",
fill = "#FFA319FF")Show the code
ggarrange(pctFunctional,pctNonFunctional,pctUnknown,
ncol = 2,
nrow = 2)adding dummy grobs
4. CORRELATION ANALYSIS
4.1 Create Data Table for Correlation Matrix Analysis
4.2 Visualise Correlation Matrix
Usage of the code chunk below :
corrplot.mixed( ) - corrplot - to use mixed methods to visualise a correlation matrix.
This plot allows to identify the pattern and the relationship in the matrix.
Show the code
corrplot.mixed((cor(cluster_vars[,2:12])),
upper = "number",
lower = "ellipse",
tl.col = "black",
diag = "l",
tl.pos = "lt")
Remarks :
Following are the pairs with strong correlation :
| correlation coefficients | variable_1 | variable_2 |
|---|---|---|
| 1.00 | pct_mechPump | pct_uc1000 |
| 0.99 | pct_tapStand | pct_uc250 |
| 0.99 | pct_uc300 | pct_ucN1000 |
| -0.91 | pct_mechPump | pct_ucN1000 |
| -0.91 | pct_uc1000 | pct_ucN1000 |
| -0.90 | pct_mechPump | pct_uc300 |
| -0.90 | pct_uc300 | pct_uc1000 |
4.2.1 Replace Row ID with “shapeName”
4.2.2 Trim High Correlation Variable and “shapeName”
5. CLUSTERING ANALYSIS
5.1 Hierarchy Clustering
There are four (4) main steps :
- compute proximity matrix.
- assign data point to a cluster.
- merge clusters based on similarity between clusters.
- update the distance matrix.
5.1.1 Standardise Data
As shown in the 4.2.3.1, there are few variables with Max. different from others. Hence, standardisation will be required prior to further analysis.
5.1.1.1 standardise with min-max method
5.1.1.3 visualise distribution of standardised clustering variable
-- functional water point
Show the code
fwp <- ggplot(data=cluster_varsTrim,
aes(x= `pct_functional`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
fwp_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_functional`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
fwp_z <- ggplot(data=fwp_zDf,
aes(x=`pct_functional`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(fwp, fwp_std, fwp_z,
ncol = 3,
nrow = 1)
Show the code
fwp <- ggplot(data=cluster_varsTrim,
aes(x= `pct_functional`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
fwp_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_functional`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
fwp_z <- ggplot(data=fwp_zDf,
aes(x=`pct_functional`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(fwp, fwp_std, fwp_z,
ncol = 3,
nrow = 1)
-- water point deployed with handpump
Show the code
HP <- ggplot(data=cluster_varsTrim,
aes(x= `pct_handPump`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
HP_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_handPump`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
HP_z <- ggplot(data=fwp_zDf,
aes(x=`pct_handPump`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(HP, HP_std, HP_z,
ncol = 3,
nrow = 1)
Show the code
HP <- ggplot(data=cluster_varsTrim,
aes(x= `pct_handPump`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
HP_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_handPump`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
HP_z <- ggplot(data=fwp_zDf,
aes(x=`pct_handPump`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(HP, HP_std, HP_z,
ncol = 3,
nrow = 1)
-- water point with 1000 users usage capacity
Show the code
uc1000 <- ggplot(data=cluster_varsTrim,
aes(x= `pct_uc1000`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
uc1000_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_uc1000`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
uc1000_z <- ggplot(data=fwp_zDf,
aes(x=`pct_uc1000`)) +
geom_histogram(bins=20,
color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(uc1000, uc1000_std, uc1000_z,
ncol = 3,
nrow = 1)
Show the code
uc1000 <- ggplot(data=cluster_varsTrim,
aes(x= `pct_uc1000`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Before Standardisation")
fwp_stdDf <- as.data.frame(nga_wpStd)
uc1000_std <- ggplot(data=fwp_stdDf,
aes(x=`pct_uc1000`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Min-Max Stdsn.")
fwp_zDf <- as.data.frame(nga_wpZ)
uc1000_z <- ggplot(data=fwp_zDf,
aes(x=`pct_uc1000`)) +
geom_density(color="black",
fill="steelblue") +
ggtitle("Z-score Stdsn.")
ggarrange(uc1000, uc1000_std, uc1000_z,
ncol = 3,
nrow = 1)
5.1.2 Compute Proximity Matrix
5.1.3 Compute Hierarchical Clustering
Usage of the code chunk below :
hclust( ) - stats - to compute cluster with agglomeration method.
ggdendrogram( ) - ggdendro - to plot dendrogram with tools available in ggplot2.
Show the code
hieClust_warD <- hclust(proxmat_euc, method = 'ward.D')
ggdendrogram(hieClust_warD,
rotate = TRUE,
size = 2,
theme_dendro = FALSE)
5.1.4 Determine Optimal Clustering Algorithm
5.1.5 Determine Optimal Clusters
To determine the optimal clusters to retain, following commons methods are tested :
Gap statistic
Elbow
Average Silhouette
5.1.5.1 compute Gap Statistic method
-- visualise gap_stat
Usage of the code chunk below :
fviz_nbclust( ) - factoextra - to compute and visualise the Optimal Number of clusters.
Show the code
set.seed(12345)
fviz_nbclust(nga_wpZ,
kmeans,
nstart = 25,
method = "gap_stat",
nboot = 50)+
labs(subtitle = "Gap statistic method")Warning: did not converge in 10 iterations
Warning: did not converge in 10 iterations
Warning: did not converge in 10 iterations
5.1.5.2 compute and visualise Elbow method
Show the code
fviz_nbclust(nga_wpZ, kmeans, method = "wss") +
geom_vline(xintercept = 4, linetype = 2)+
labs(subtitle = "Elbow method")
5.1.5.3 compute and visualise Silhouette method
Show the code
fviz_nbclust(nga_wpZ, kmeans, method = "silhouette")+
labs(subtitle = "Silhouette method")
Remarks :
Given the Elbow method, Silhouette method and Gap Statistic method, the 5-cluster by Silhouette method will be used for the rest of the study.
5.1.5.4 interpret with Dendrogram
Usage of the code chunk below :
rect.hclust( ) - stats - to draw the dendrogram with a border around the selected clusters.
Show the code
plot(hieClust_warD, cex = 0.6)
rect.hclust(hieClust_warD,
k = 5,
border = 2:5)
5.1.6 Visually-Driven Hierarchical Clustering Analysis
The data is loaded into a data frame, but it has to be a data matrix to plot the heatmap. Hence, the data frame will need to first transform into a matrix.
5.1.6.1 transform data frame into matrix
5.1.6.2 plot interactive cluster heatmap
5.1.6.3 map the formed cluster
5.1.6.4 append groups to wp_ngaTrans
5.1.6.5 plot choropleth map :: nga_clust.sf
Show the code
qtm(nga_clust.sf, "cluster")
Remarks :
The choropleth map above shows the fragmented clusters by the used of non-spatial clustering algorithm (hierarchical cluster analysis method).
5.2 Spatially Constrained Clustering :: SKATER Approach
SKATER function only support sp objects in SpatialPolygonDataFrame. Hence, the wp_ngaTrans has to first transform into SpatialPolygonDataFrame before proceed further.
5.2.1 Convert SF to SP Data Frame
5.2.2 Compute Neighbour List
5.2.2.1 remove 0-neighbour region
5.2.2.2 plot Neighbour List by Centroid Node
Usage of the code chunk below : plot the boundary first before the neighbour list object to avoid any region from being clipped away.
Show the code
plot(nga.sp1,
border=grey(.5))
plot(nga.nb1,
coordinates(nga.sp1),
col="blue",
add=TRUE)
5.2.3 Compute Minimum Spanning Tree (MST)
5.2.3.1 calculate edge costs
5.2.3.2 specify spatial weight
5.2.3.3 compute minimum spanning tree
5.2.3.4 plot MST Neighbour List
Show the code
plot(nga.sp1, border=gray(.5))
plot.mst(nga_minSpanT,
coordinates(nga.sp1),
col="blue",
cex.lab=0.7,
cex.circles=0.005,
add=TRUE)
5.2.4 Compute Spatially Constrained Cluster
5.2.4.1 tabulate cluster assignment
5.2.4.2 plot the pruned tree
Show the code
plot(nga.sp1, border=gray(.5))
plot(clust5,
coordinates(nga.sp1),
cex.lab=.7,
groups.colors=c("red","green","blue", "brown", "pink"),
cex.circles=0.005,
add=TRUE)Warning in segments(coords[id1, 1], coords[id1, 2], coords[id2, 1],
coords[id2, : "add" is not a graphical parameter
Warning in segments(coords[id1, 1], coords[id1, 2], coords[id2, 1],
coords[id2, : "add" is not a graphical parameter
Warning in segments(coords[id1, 1], coords[id1, 2], coords[id2, 1],
coords[id2, : "add" is not a graphical parameter
Warning in segments(coords[id1, 1], coords[id1, 2], coords[id2, 1],
coords[id2, : "add" is not a graphical parameter
Warning in segments(coords[id1, 1], coords[id1, 2], coords[id2, 1],
coords[id2, : "add" is not a graphical parameter
5.2.5 Visualise SKATER Clusters in Choropleth Map
Show the code
groups_mat <- as.matrix(clust5$groups)
nga_spClust.sf <- cbind(nga_clust.sf1, as.factor(groups_mat)) %>%
rename(`sp_cluster`=`as.factor.groups_mat.`)To compare the output of hierarchical clustering and spatially constrained hierarchical clustering :
Show the code
hieClust_map <- qtm(nga_clust.sf1,
"cluster") +
tm_borders(alpha = 0.5)
ngaClust_map <- qtm(nga_spClust.sf,
"sp_cluster") +
tm_borders(alpha = 0.5)
tmap_arrange(hieClust_map, ngaClust_map,
asp=NA, ncol=2)Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).
Warning: One tm layer group has duplicated layer types, which are omitted. To
draw multiple layers of the same type, use multiple layer groups (i.e. specify
tm_shape prior to each of them).
5.3 Spatially Constrained Clustering :: ClustGeo Method
5.3.1 Perform Ward-like Hierarchical Clustering
Usage of the code chunk below :
hclustgeo( ) - ClustGeo - to perform a typical Ward-like hierarchical clustering.
Show the code
proxmat_ngc <- dist(cluster_varsTrim1, method = 'euclidean')Show the code
nonGeo_clust <- hclustgeo(proxmat_ngc)
plot(nonGeo_clust, cex = 0.5)
rect.hclust(nonGeo_clust,
k = 5,
border = 2:5)
5.3.1.1 visualise the formed clusters
Show the code
groups_ngc <- as.factor(cutree(nonGeo_clust, k=5))Show the code
nga_ngeo_clust.sf <- cbind(wp_ngaTrans1, as.matrix(groups_ngc)) %>%
rename(`cluster` = `as.matrix.groups_ngc.`)Show the code
qtm(nga_ngeo_clust.sf, "cluster")
5.3.2 Perform Spatially Constrained Hierarchical Clustering
5.3.2.1 determine alpha value
choicealpha( ) - psych - to determine a suitable value for the mixing parameter alpha.
Show the code
cr <- choicealpha(
proxmat_ngc,
dist_mat,
range.alpha = seq(0, 1, 0.1),
K=5,
graph = TRUE)
Remarks :
With reference to the plot above, alpha = 0.4 to be used to perform spatially constrained hierarchical clustering.
5.3.2.2 compute spatially constrained hierarchical clustering
5.3.2.3 derive cluster object
5.3.2.4 combine group_cg with wp_ngaTrans1
5.3.2.5 plot delineated spatially constrained cluster
Show the code
qtm(wp_nga1_GClust, "cluster")
5.4 Visual Interpretation of Clusters
5.4.1 Visualise Individual Clustering Variable
5.4.1.1 plot boxplot
Show the code
ggplot(data = nga_ngeo_clust.sf,
aes(x = cluster, y = pct_functional)) +
geom_boxplot()
Remarks :
The boxplot reveals Cluster 5 displays the highest mean of functional water points. This is followed by Cluster 1, 3, 2, and 4.
5.4.2 Visualise Multivariate
Show the code
ggparcoord(data = nga_ngeo_clust.sf,
columns = c(2:19),
scale = "globalminmax",
alphaLines = 0.2,
boxplot = TRUE,
title = "Multiple Parallel Coordinates Plots of Variables by Cluster") +
facet_grid(~ cluster, scales = "fixed") +
theme(axis.text.x = element_text(angle = 30))
The parallel coordinate plot above reveals that households in Cluster 4 townships tend to own the highest number of TV and mobile-phone. On the other hand, households in Cluster 5 tends to own the lowest of all the five ICT.
Note that the scale argument of ggparcoor() provide several methods to scale the clustering variables. They are:
std: univariately, subtract mean and divide by standard deviation.
robust: univariately, subtract median and divide by median absolute deviation.
uniminmax: univariately, scale so the minimum of the variable is zero, and the maximum is one.
globalminmax: no scaling is done; the range of the graphs is defined by the global minimum and the global maximum.
center: use uniminmax to standardize vertical height, then center each variable at a value specified by the scaleSummary param.
centerObs: use uniminmax to standardize vertical height, then center each variable at the value of the observation specified by the centerObsID param